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1 . 


INTRODUCTION 


This report consists of a technical explanation of the code 
used to calculate the steady supersonic, three-dimensional, 
inviscid flow over blunt bodies. Because this code has been 
extensively documented elsewhere (e.g., refs. 1 and 2), a cer- 
tain measure of brevity is possible; further details are available 
in these references. 

In the following section (section 2) , a brief discussion 
of the theoretical and numerical formulation of the problem is 
given, including exposition of the boundary and initial conditions 
used. In section 3, the overall flow logic of the computer pro- 
gram is given, in section 4, the program usage and operation are 
described, and in section 5, the program accuracy and limitations 
are discussed. 


2. THEORETICAL AND NUMERICAL FORMULATION OF THE PROBLEM 
2.1 Review of the Governing Equations 

The fluid dynamic equations in conservation-law form govern- 
ing steady, inviscid, three-dimensional, compressible flow of a 
nonheat-conducting gas can be written in cylindrical coordinates 
as follows; 


3U/3Z + 3F/3r + 


3G/3(t) + H = 0 


( 1 ) 
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where jj, F, G, and H are four-component vectors defined as 
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Here p and p represent dimensional pressure and density and u, v, 
and w denote velocity components in the coordinate directions z, r, 
and (p. The nonlinear system (1) of four equations represents 
conservation of mass and mcraentun. 

The governing set of equations is made complete by the addition of 
energy conservation as given by the equation for total enthalpy 

= h(p,p) + q^/2 = const (2) 

where q is the magnitude of the velocity vector and h(p,p) is the 
state equation for static enthalpy. The specific formulation for 
h depends, in particular, on whether the gas is assumed to be 
perfect or everywhere in local thermodynamic equilibrium. Explicit 
representations for h are described later. 

The vehicle geometry and the location of the outer or peripheral 
shock surface are represented by functions of the form 

Tb = rb(z,<i)) , = r^(z,6) (3) 
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The function r, is specified (sec. 2.4) and r is determined during 
the course of the numerical computation. As is common practice in 
problems of this type, the distance between the body and peripheral 
shock is normalized by a transformation of the radial variable r. 
This yields a rectangular computational plane whose boundaries 
consist of the plane of symmetry and the body and shock surfaces 
as shown in the sketch below. 



I 

$ , n=0° 

I 

z = constant 
Discretized Flow Region 



Computational Plane 


Mesh Description 


Since the flow variables can vary rapidly in the cross flow plane, 
an independent variable transformation is performed in the !p direc- 
tion to cluster points in that region of suspected large gradients. 

The transformation, which is also given below, was introduced by 
Woods (ref. 3) and has been successfully applied by Schiff (ref. 4). 


The equations of the independent variable transformations are 

z = z, ^(z,r,(p) = (r - ~ ^ ~ ^(xtan 0) (4) 

where k is a free parameter with the range 0 ^ k i 1. The points 
are clustered about the wing tip region (90° plane) for small values 
of K and the spacing approaches uniformity as k approaches unity. 
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The change of variable procedure outlined above is applied to each 
vector of Eq . (1), and the resulting terms are then rearranged in 
conservation-law form to yield the following equation: 


9U/9Z + 9F/95 + 9G/9n + H = 0 


(5) 


where 

U 

{-irfa + - r )1U + # - [r + 5 (r - r )lG)/(r^ - 

z z z 4> 4> (p 

2 2 2 ~ 

(k cos n + sin n)G/K 

I + >2 + (r - r )6)/(r^ - r^) 

Z Z ({) (p 

- ( 1 - ) sin ( 2r\ ) G/k 

The finite-difference analogue of Eq. (5) is integrated with 
respect to the hyperbolic coordinate z to yield values of the con- 
servative variable JJ. Subsequent to each integration step, the 
physical flow variables p, p, u, v, and w must be decoded from the 
components u^ of U, This necessitates the solution of five simul- 
taneous nonlinear equations consisting of Eq . (2) together with the 

four elements . The velocity components v and w are easily found 
and are given by 


U = 

F = 

G = 
H = 


V = u^/u^, w = u^/u^ 


( 6 ) 


If the along with Eqs . (6) are used to eliminate the explicit 

dependence of p, p, v, and w from Eq . (2), one obtains the following 

implicit expression for the velocity, u: 
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(7) 


E(u) = u^/2 + h[p(u),p(u)] - r/2 = 0 

where 

2 2 2 

p(u) = U 2 - u^u, p (u) = u^/u, and r = 2H^ ~ 

The decoding procedure is now reduced to a problem of root 
finding, i.e., the z-velocity component u that satisfies Eq . (7). 

Two roots exist; one corresponds to subsonic flow and is discarded 
since in the results presented u is always supersonic, and the other 
corresponds to supersonic flow and gives the desired solution. The 
procedure for solving Eq . (7) depends on whether a perfect or real 

gas is being considered and consequently on the function h(p,p) . 

For a perfect gas h(p,p) is simply related to pressure and density 
and when combined with Eqs. (7 and 8) yields a quadratic equation 
that can be solved to find an analytical representation (ref. 5) 
for the supersonic velocity u. 

For a real gas no such simple explicit functional relationship 
exists. The conventional procedure (refs. 6 and 7) for evaluating 
real gas state relations is to use a combination table lookup and 
curve- f i tt i ng procedure. Such a scheme is adopted here. A 
particularly rapid Fortran-language computer code called RGAS is 
available that returns values of static enthalpy h, speed of sound, 
temperature, and entropy, with either pressure and density or 
pressure and entropy as independent variables. 

To find the roots of Eq . (7) for a real gas, a root finding 

algorithm is employed. The "successive linear interpolation" scheme 
described by Dekker (ref. 8) was found to be particularly efficient. 
Slightly more than seven iterations are required, on the average, 
to find the desired supersonic root for a given U (three iterations 
is the happenstance absolute minimum). 
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Three-dimensional numerical integration methods (see, e.g., 
refs. 6, 7, and 9) which employ physical rather than conservative 
dependent variables are generally able to include state relations in 
a more direct manner since a decoding procedure is not necessary. 
Only half again more computational time than a perfect-gas calcula- 
tion is required for those codes compared to a factor of about four 
when conservative variables are employed. However, the ability to 
capture shocks by use of conservative variables is worth the addi- 
tional cost. 


2.2 Boundary Conditions 

In this section the boundary condition schemes applied at the body 
and shock surfaces are discussed briefly. At the body the surface 
tangency condition 


2 • lib = ° 


(9) 


is imposed v/here q is the velocity vector 


a = ux + vx + wx, 

— — z — r — <p 

and n is the outward unit normal to the body. The second-order 
— b 

predictor-corrector scheme, MacCormick's predictor and Eq . (15) 

below, is first applied at the body to yield the conservative 

variables. These variables are then decoded [see discussion 

following Eq. (6)] to obtain p-, ^ , Pn ^ , v-, ■, and w, . at z 

r J f J ' -*-/J r J 

In general the resulting velocity vector based on the predicted 

velocity components will not satisfy the surface tangency condition 
Eq. (9) and, in fact, will be rotated out of the surface tangent 
plane by a small angle A0 . In applying Abbett's method (ref. 10) 
this angle can be determined from the following equation: 


■ -1 f-n+1 

A9 = sxn j 


_n+l'| 


( 10 ) 
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where is defined to be the magnitude of 

If: -1 , : 

unit normal n^ to the body can be calculated from 


and the outward 
(see ref. 5) 





|-r2^+ 1 + 


( 11 ) 


The function f, = r, - r, (z,*) = 0 describes the body surface 
b b b ^ 

[Eq. (3)]. Once A4> is determined the velocity vector is rotated 
back into the surface tangent plane by imposing a simple compression 
or expansion wave. 


If A(|) is positive, then an expansion is necessary for the rotation 
„ 1 1 

of ^ . and if A0 is negative a compression wave is required. The 

1 f : 

corrected value of the static pressure is found from the integral 
relation (ref. 12) for the Pr andtl-Meyer turning angle v(p;H ,s) 
which depends on pressure and has the total enthalpy and entropy 
as parameters. The corrected value of pressure is found by solving 


V (p 


n+1 
If j' 


s ) = 


n+1 


H^,s) 


+ A0 


( 12 ) 


for the pressure p, . . In this equation AS is given by Eq . (10). 

f : 

If A0 is sufficiently small, Eq . (12) can be inverted and solved 

analytically for p-, • only in the case of a perfect gas (ref. 5) . 

For a real gas, Eq . (12) can be inverted by the use of a table- 

lookup method. The isentropic flow assumption underlying Abbett's 
boundary condition procedure requires that the table be generated 
only once at the very beginning of a flow field calculation when the 
entropy on the body stream surface is known. The table elements are 
pressure and Pr and tl-Meye r turning angle v. The procedure for 
generating the table is described by Hayes and Probstein (ref. 12). 
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.n+1 


The pressure p, . appearing in Eq. (12) is determined by first 
—n+1 

finding v(p, H , s) from the table with the predicted pressure 
—n+1 ^ ^ 

P-, . as the argument. The angle A0 given by Eq . (10) is then added 

+ f J 


n+1 


to the result and the desired value for the corrected pressure p'^^ , 
is then found from the same table with H , s) + A0 as the 


,n+l 


, n+1 


argument. The remaining flow variables y 
are then obtained in the same fashion as described in reference 5. 


V- . f and V7_ 

1.: i,D 


The outer boundary of the computational plane is the peripheral 
shock wave which completely encompasses the disturbed flow region. 

In this report a procedure is used whereby the peripheral shock 
wave (i.e., the shock for which freestream conditions are maintained 
on the upstream side) is treated as a discontinuity. The Rankine- 
Hugoniot relations are satisfied exactly across this discontinuity. 

There are many ways of incorporating a sharp-shock calculation in a 
numerical algorithm, including method of character istics and simpler 
methods such as those of Kentzer (ref. 11), Barnwell (ref. 13), 
and Thomas, et al. (ref. 14). Thomas' scheme is used here since it 
is relatively simple to implement and, furthermore, comparisons with 
other methods, including a full method of cha r ac te r i s t i cs , have 
shown it to yield sufficiently accurate results for our purpose. 

The pressure downstream of the shock wave is the basic variable on 

which all other shock-wave variables depend. Its estimated value 

2 < j < N at (n + l)Az, is first found by the predictor step 

[Eq. (16a) below]. Since all other variables associated with the shock 

wave can be expressed as functions of pressure through the Rankine- 

Hugoniot equations, their values can then be found. The pressure 
X 

Pj^j- . is then recomputed by the corrector, (Eq. (16b) below), and the 

S/ J 

dependent flow variables are adjusted accordingly. 

The boundary conditions at the planes of symmetry are applied in 
the conventional manner whereby the fringe points (see previous sketch) 
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are filled with data reflected across the planes of symmetry, e.g.. 


Pi,l P±,3' 


w . 

1,1 




e tc . 


2.3 Initial Conditions 


It is necessary to have an axis-normal starting plane to begin 
the calculation. All variables, the shock location and the shock 
slopes must be specified on this plane. There are two ways of 
obtaining this starting solution. The first is using the pointed- 
cone starting solution procedure which is contained internally within 
this code. The second is for blunt bodies; here, the starting-plane 
solution is obtained from a separate computer code such as Barnwell 
(ref. 13) or Moretti (ref. 15). 


2.4 Finite-Difference Scheme 


Since only the peripheral shock is treated as a sharp dis- 
continuity and the others are "captured" by the difference 
algorithm, of prime importance is the selection of the finite- 
difference scheme to be used to advance the field points, i.e., the 

points for which 2 < i < N_ - 1 and 2 < j < N in the sketch above. 

= =5 =-’=n 

It has been found (ref. 17) that MacCormack's scheme (ref. 16) , 
which is an accurate predictor-corrector scheme, is the most 
efficient second-order algorithm to use in a shock-captur ing 
techn i que . 

The algorithm can be written as 
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( 14 ) 
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Az 
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I — 1 

e> 

- J 

- Azh[^[ 

1-1 '3 

H 

1 

H 

1 

-i/3j 


. = U(nAz,iAC, jAn) 

1 / J 

. = F(U^ .,nAz,iAC,jAn) 

^ r J “"i / J 

F^^^ = Ffuf^J, (n+l)Az,iAC,jAn 
-*- / J L / J 


At the body (i = 1,2 ^ j <_ N^) Abbett's fref . 10) scheme is 
used to satisfy the surface tangency condition as is discussed 
above. It relies on information provided by the finite-difference 
scheme. The numerical algorithm used for the field points, hov/ever, 
cannot be used on the surface since it requires points on either 
side of the point being advanced and thus data at a set of points 
that would lie within the body. Consequently, a special second-order 
accurate algorithm was constructed which requires data only on or 
outside the body. This scheme uses the predictor step of -^acCo rmack ' s 
method (ref. 16), Eq . (13) followed by the corrector step given by 


U . . 
-1/3 


1 

2 


.• 

■1/3 


-1/3 


^ f (1) _ (1)' 

AC [-1+1, j -i/3, 


An [-1,3 -x,]-lj 


- AzH 


( 1 ) 
= 1/ j 


Az 

AC 


[-1+2, j 


.n 


- 2F . . . 

-1+1 / 3 


+ 


F^ 

- 1 , 


(15) 


where i = 1. After Ma cCo rm.ack ' s predictor and Eq . (15) have 

been used to advance the data at the body, then Abbett's scheme 
is used as a final corrector. 
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At the shock wave, a predictor-corrector sequence is again 
used and requires data at the shock and one point below it. The 
algorithm is as follows; 


( 1 ) 
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Az 
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- G 
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[-1/3 

-1/3-iJ 

-1/ 3j 


(16b) 


where i = N^. Equations (16a) and (16b) are used in conjunction with 
the Rankine-Hugoniot relations as described above, to determine the 
peripheral shock slope. 


The integration stepsize is based on the analysis presented in 
reference 5, in which formulas for the projections of the slopes of 
the characteristics in the z-^ and z-<t> planes are given and, with a 
slight modification due to the clustering transformation of 4> , can be 
used to determine the maximum possible stepsize. 


3. OVERALL PROGRAM FLOW LOGIC 


The overall logical flow of the program is described in this 
section. PROGRAM MAIN is the executive program for this computer 
code. It controls the main flow of the program logic and it starts 
by calling SUB.GE0M3 which reads the geometry description. Next 
SUB. INPUT is called from which the flow conditions and the operating 
controls are read. Program MAIN then calls SUB.IMITA. This 
initializes the starting plane and defines various physical con- 
stants. To complete the ini tial ization process SUB.3NDRY(2) is 
called to fill in the date at the planes of reflection. This 
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completes the process of getting the code ready to march downstream 
from an axis normal initial data plane. 

SUB. EIGEN is called to calculate the initial marching stepsize 
based on the grid size and the values of the flow variables. The 
program then prints out the entire flow field and outer shock 
structure. This is done by calling SUB. OUTPUT. 

The marching process is done by looping from 1 to NITER and 
is controlled by calling two subroutines, EIGEN and DIFFR. The 
stepsize of the marching process is determined every ICQNST(49) 
iterations by calling SUB. EIGEN. Once the stepsize is determined 
SUB. DIFFR is called to perform MacCormack's pred iction-correction 
sequence. Within SUB. DIFFR flow field variables on the boundaries 
are calculated by SUE.BNDRY and SUB. SHOCK. 

After completing the iteration loop, SUB. OUTPUT is called to 
furnish a final solution. 

4. PROGRAM USAGE AND OPERATION 

The input data options for the code are described below in terms 
of the card number, its format, and the variables v;hich are defined 
by that card. All geometry inside the program is expressed in terms 
of body-oriented coordinates, z (longitudinal), r (radial, normal to 
z) , and 4> (radial angle around the z-axis counterclockwise, starting 
at zero radians straight down), or functions and transforms thereof. 
The unit for z, r and all lengths input are arbitrary, a user's 
choice, requiring only that they be consistent throughout each case. 
In the case of angular dimensions, it has been found more convenient 
to input them in degrees and let the program convert and store them 
in radians. 
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Card 
No . 


Format 


Variables 


1 

815 

NSEG 

# of Segment Points 



KIND 

Flag for Kind of Segment 

0 = Sphere or Circulr Ogive 

1 = circular Cone 

2 = Circular Cone with Flat 

Cut 

2 

8F10.6 

ZSEG 

Z - Station Initiating Seg- 
ment 

3 

8F10.6 

RSEG 

r - Coord. Initiating Segment 

4 

8F10.6 

DSEG 

Distance from Centerline 
to Flat Chord, Initiating 
Segment 

5 

8F10.6 

ASEG 

Angle Between Straight Down 
and DSEG Initiating Segment 

5a 

3F10. 6 

ZC 

RC 

RADIUS 

Z at Center of Longitudinal 
Arc 

r at Center of Longitudinal 
Arc 

Radius of Longitudinal Arc 


Card 1 (Format 815) specifies NSEG (=number of segment points) in 
its first field, and KIND flags in the next NSEG fields (up to 7 
maximum). The segment points, whose dimensional specifications are 
punched on cards 2 through 5, bracket NSEG-1 segments. The last 
point usually initiates an infinite extension of the last bracketed 
segment . 

The KIND flags tell v;hat kind of contour each segment point 
initiates, where 0 = Circular Arc with Circular Cross Section, 

1 = Straight Line with Circular Cross Section (Forecone, Cylinder, 
or Aftcone) , and 2 = Cut Cone (Circular Cross Section Truncated 
by a Flat Cut). It should be noted that KIND = 0 is not restricted 
to spheres, but can also be used for circular arc ogives of circular 
cross-section. It has been used, for example, in rounding off 
contours approaching a planar axis-normal base. In such a case 
the KIND = 0 flag initiating the last segment is also the last 0 
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flag, to avoid calling an extra card later. A sphere-cone-arc 
would have NSEG = 4 and KIND = 0,1, 0,1. 

The next four cards contain data describing successive cross 
sections (up to 7) in successive fields of each card, using Format 
7F10 (with the punched decimal points governing the decimal point 
location in each field) . With the body segmented into up to six 
physical contours of types KIND, each bounded by a pair of E-planes, 
the z-plane coordinates called ZSEG are specified on card 2. The 
radius of each cross-section, RSEG, is specified in corresponding 
card fields on card 3. If a planar surface is to truncate any 
circular cross section, the radial distance from the centerline 
to the midpoint of the chord formed is DSEG on card 4, for those 
fields where the planar surface is on the vehicle. 

The orientation of a planar surface cut is not limited to the case 
of d = 0. Card 5 provides specification of ASEG, the angle p in 
degrees subtended by DSEG in each z-plane of appl icabi 1 i ty . To 
define a plane, ASEG should be constant. 

Card 5 contains circular-arc data, and is read only if a KIND = 0 
flag appears within the fields activated by NSEG on card 1. More 
than one longitudinal circular arc may be specified, and one card 
must be read for each KIND = 0 segment initiated on card 1. The 
three data items of card 5a are 2C (Longi tud inal Loca t ion of 
Circular-Arc Center) , RC (Radial Location of Circular-Arc Center) , 
and RADIUS (Long i tud inal -C i r cul a r-Arc Radius). It is important 
that the data of cards 2 through 5a be accurate to at least 5 
significant figures in order that successive contours meet smoothly 
at points of tangency such as sphere-cone or cone-arc junctures. 

It should be noted that geometry spec i f ica t i ons usually should 
encompass all regions of the configuration bracketed by the z value 
at the input plane and the final ZMAX specified on a later control 
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card. However, the last ZSEG specified may be less than ZMAX if the 
kind of contour beyond the last ZSEG continues as a qual i tat ively 
unchanged extension of the NSEG-1 contour bracketed by the last ZSEG. 

On the other hand, if the input plane z-station is on a cone downstream 
of a sphere-cone juncture or is exactly at such a juncture, the first 
geometric contour described may be the cone (if the user so desires) 
and not the sphere, and that cone may be described by use of its apex 
rather than by the sphere-cone juncture as the first ZSEG. Finally, 
if a configuration is so variegated as to consist of more than six 
different segments (such as sphere-cone-cylinder-flare-cone-cylinder- 
aftcone-arc) , such a case can be run by punching out a restart plane 
near the end of the sixth segment and starting a new case with the 
sixth segment geometry respecified as the first of the aft-region 
series. 

Only eight basic cards are needed to input ordinary case data and 
controls. There are several options which also exist that either 
require additional cards or not as many cards. The ordinary cards 
and some options will be described following the glossary of the 
input cards. 


Card 

No ■ Format 

(Cards 6-12 are read in SUB, INPUT) 

6 3E15.6,5X,I5 XMACH 

ALPHA 

GAMMA 

NREAL 

7 3F10.5 PHIFD 

RK 

RJ 

3 515 NIT 

NIPHI 

NITER 


Variables 


Mach number 

angle of attacli (degrees) 
ratio of specific heats 
0 for perfect gas, -1 for 
real gas (pointed cone 
starting solutions are 
generated internally for 
perfect gas option only) 

meridional angle about 
which points are clustered 
meridional clustering para- 
meter (0 for no clustering) 
radial clustering parameter 
(0 for no clustering) 

No. of points between body 
and shoc!< (max = 20) 

No. of intervals in meri- 
dional direction (max = 36) 
No. of integration steps 
desired (when ZEND is 
specified set NITER to 99999) 
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ICONST(49) 


9 


10 


11 


3F10.5 


515 


2F10.5, 315 


NCONE 

CONST (9) 
CONST (4) 
CONST ( 5 ) 


DISK 1 


DISK 2 


TAPE 1 


TAPE 2 


NTDSOS 

ZBS 

ZFLD 

ITPRTB 


ITPRTF 

NCASE 


Stepsize is computed every 
ICONST(49) iterations (5 is 
nominal) 

1 for pointed cone solutions, 

2 for all other geometries 

Courant No. (usually 0.9) 
Radial dissipation constant 

Meridional dissipation con- 
stant 

1 reads solution from tape, 

2 writes solution on tape, 

3 does nothing (logical 
unit 12) 

1 reads solution from tape, 

2 writes solution on tape, 

3 does nothing (logical 
unit 11) 

1 does nothing, 2 stores 
body shape and writes data 
on tape each Z station, 3 
writes data only (logical 
unit 9) 

1 does nothing, 2 reads 
starting solution from 
punched cards, 3 stores 
solution on punched cards 
when exiting (logical unit 
7). If TAPE2 = 1 and DISK 1 
and DISK 2 = 2 or 3 a pointed 
cone solution will be gen- 
erated for the perfect gas 
case only 

0 

increment in z for ') 

printing shoc)c and | 

body variables (ZBS I 

> ZEND if not desired) > Print 
increment in z for 
printing field vari- 
ables (ZFLD > ZEND if 
not desired) 

No. of iterations for 
printing shocl< and 
body variables 
(ITPRTB > NITER if not 
desired) 

No. of iterations for 
printing field vari- 
ables (ITPRTF > NITER 
if not desired) 

If > 0, new case follows 


Based on 
Z Stations 


Print 
\ Based on 
I Number of 
I Iterations 



(The following card contains values used in force and moment calcula- 
tions or in shifting the origin of the pointed cone starting solution.) 


12 5F10.5,I5 


DIAM 

ALENGT 

ZREF 

ZCG 

ZSHIFT 

IFANDM 


length used in calculating refer- 
ence area; usually maximum 
diameter 

reference length used in calculat- 
ing moments 

moment reference center 
center of gravity location for 
static margin calculation 
the value of Z which corresponds 
to the starting cone origin; if 
no shift set = 0 

0, force and moment calculation 

1, no force and moment calculation 
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(If starting solution is to be read from punched cards (TAPE 2 = 2 ), 
the following three cards are read in main program. If solution is 
read from magnetic storage device, these are not required.) 

12a 5E15.6 X.MACH, ALPHA, GA.MMA, RK, PHIFD (defined above) 

12b 5E15.6 RJ (defined above) 

12c 3I5,4E15.6 NIT, NIPHI, NREAL, (defined above) 

PlINF free stream pressure, > Real Gas 

RIINF free stream density. Option 

VlINF free stream velocity, I Only 

GASCON gas constant (1716.0 (Dimen- 

for air) J sional) 

[If NREAL = -1, gas tables are place here and will be read in SUB.RGAS 
(523 cards for equilibrium air)] 

(If TAPE 2=2 punch card starting solution is placed here. The first 
card is the Z station of the starting plane and is followed by flow 
variables at each node.) 

Card 6 (FORMAT 3E 1 5 . 6 , 5X , I 5 ) contains the basic free-stream flow 
velocity, angle of attack, gamma, and the operating control for the 
type of gas being employed. For ideal gas cases, (ratio of 

specific heats) is held constant throughout the flow field. The free 
stream Mach number, XMACH, is also a key parameter in the calculation 
of dimensionless properties throughout the flow field. It should be 
noted that the angle of attack, ALPHA, is to be input in degrees. 
NREAL is the flag specifying the gas type, 0 for ideal gas and -1 
for equilibrium real gas. Pointed-cone starting solutions are gen- 
erated for the ideal gas case (NREAL=0) only. 

Cards 7 through 12 contain the operating controls of the 
program. Card 7 determines the meridional location and the amount 
of clustering both meridionally and radially. The default values 
are zero(O) . Card 8 contains the grid size, NIT and NIPHI, the 
maximum number of iterations, NITER, the increment between stepsize 
determinations ICONST(49), and the body-surface boundary condition, 
NCONE. COMST(9) is a Co ur an t-numbe r factor governing the ratio of 
stepsize actually used to that allowable by stability criteria and 
it is found on Card 9. Also, the fourth-order smoothing parameters, 
COMST(4) , and CCMST(5) are defined on this card. If needed, these 
parameters are typicaly on the order of 0.1. Card 10 contains the 
parameters that are used for storing or reading a solution from a 
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peripheral device such as a tape or punched cards. Card 11 
determines the print controls for this program and NCASE which is 
set to a value greater than zero if cases are to be run in 
succession. Card 12 defines various force and moment parameters and 
the value of z which corresponds to the starting-cone origin. The 
default for ZSHIFT is 0.0. 


Cards 12a, b, and c are read only if a storting solution is read 
from punched cards (TAPE=2 on card 10) ; these are read in the main 
program . 

The gas tables for equilibrium air are read in SUB.RGAS if 
NREAL = -1 in this location, followed by a punched card starting 
solution if TAPE2=2. The first card is the z-station of the 
starting plane. 


Next are cards 13 and 14 which are used to change the preceding 
program control variables at a specified z station and to initialize 
the forces and moments. If no altering is desired, one card is 
required and the z-station should be set to a value which is greater 
than ZEND. Card 14 is necessary only if IFANDM = 0 and NCOME = 2 and 
precedes the first ZALTER card (Card 13). Any other ZALTER cards 
(Card 13) are placed after card 14. 

Card 

No ■ Format Variables 

(The following card(s) is used to change the program control 
variables at preselected longitudinal (Z) stations and is read in 
Program MAIN. At least one card is required if no modifications 
are asked for. In this case ZALTER should then be > ZEND.) 


13 


F10.5,I2,I3, 

6F10.5,I2,I3 


ZALTER 

NITA 

NIPHEA 

RJA 

RKA 

PHFDA 

STP 


DISSl 

DISS2 

NSWCHl 

NSWCH5 


Z station where altering occurs 

new NIT 

new NIPHI 

new RJ 

new RK 

new PHFD 

0, stepsize determined automatically 

>0, value of desired constant step- 
. size 

new CONST (4 ) 
new CONST (5) 

0, new MacCormack 

1, old MacCormack 

0, no entropy relaxation 
,1, entropy relaxation 
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(The following card is used to initialize the force and moment calcu- 
lations, and is read in SUB.COMPOT. This card is needed only if 
IFANDM = 0. *If NCONE = 2 and IFANDM = 0 this card is read before the 
first card 13 otherwise it is read after all card 13' s.) 


14 6F12.8 


FTX 

FTY 

FTZ 


initial plane forces in the z, r, 
p direction 


RMTX 

RMTY 

RMTZ 


initial plane moments in the z, r, 
p direction 


There are error messages which are printed to aid the user when 
problems occur. The messages, locations, and their explanations 
follow below. 


SUB , BNDRY 

In this routine, there are three error checks which are denoted 
by ICHECK. If ICHECK = 1, the pressure is negative, ICHECK = 2 
the density is negative and if ICHECK = 3, the local Mach number 
squared minus 1 is less than zero [ (M^ - 1) < 0)]. These errors are 

Xj 

sometimes due to erroneous body shape. 

For each ICHECK a sequence of two error messages are printed 
out. One before the correction and one after. These are as follows: 


ERROR CHECK-NEGATIVE PRESSURE IN BNDRY 

K = Z = 

P = RHO = U = y = W = ICHECK = 

MODIFICATION INSTITUTED 

P = RHO = U = V = W = 

SUB.PMYTURN 

In this routine, if the Pr and tl-Meye r turning angle is too 
great for the flow conditions, the variables are set to values which 
correspond to a final turning Mach number of ICO. This usually 
occurs due to a rapid change of the body shape. 
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BODY TURN STOPPED AT M2 = 100.0 

SUB. EIGEN 

There are three error messages in this routine. The first is for 
a negative sound speed v/hich can occur for either a negative pressure 
or density value. The second and third messages usually occur in the 
eigenvalue calculations. These messages can occur due to the geometry 
of the vehicle or for a local axial Mach number less than one. 


ERROR CHECK -SPEED OF SOUND IN EIGEN. J = K = 

ERROR CHECK -SIGMA-EAR-1 IN EIGEN. J = ___ K = 

ERROR CHECK -SIGMA-BAR-2 IN EIGEN. J = K = 

This concludes the explanation of the major error messages. 

5. PROGPA.M ACCURACY AND LIMITATIONS 

In general, this program can be used for inviscid, supersonic/ 
hypersonic flow past bodies with a bisymmetry plane (no yav/) at 
angles of attack. Conservative estimates of applicable ranges of 
Mach number and angle of attack are M^ 2; 2 and ct _< 25. 

When using the geometry subroutine that is included, one is 
restricted to non-winged bodies. This can easily be rectified by 
changing the geometry package (SUB.GE0M3). The necessary informa- 
tion for the computer code is the body radius as a function of axial 
(z) location and meridional (c{>) location. This is furnished to the 
computer code as the radius of the body (RE) , the axial slope of the 
body radius (Rg ) and the meridional slope of the body radius (Rg ). 
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Numerically, the accuracy is second-order in time and space but 
to determine the actual physical accuracy it is best to compare with 
experimental results. These comparisons then guide the user in his 
descision as to how many grid points are needed to satisfactor ily 
resolve the flow field. Note that if the 'first radial mesh spacing 
is a fairly large proportion of the body radius, erroneous solutions 
can occur. 

The main limitations of the current version of the program are: 

1. The flow in axial direction must be supersonic (m > 1) . 

2. Yaw is not allowed. 

3. The present geometry package does not allow solution of a 
bod y wi th wings. 

4- If the Mach number, , is below 2 and at the same time the 

angle of attack is above 15°, erroneous answers can occur due 
to large radial spacings of the grid near the body. 

5. The body radius cannot be multivalued. 

6. The program treats inviscid flow only. 

The use of the program is illustrated in reference 18 by explain- 
ing the input and output for a number of test cases. 
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